home *** CD-ROM | disk | FTP | other *** search
/ Libris Britannia 4 / science library(b).zip / science library(b) / PROGRAMM / PASCAL / 0514.ZIP / CRAYZ15.ARC / VSGEFA.PAS < prev    next >
Pascal/Delphi Source File  |  1986-08-01  |  6KB  |  122 lines
  1. { Copyright (C) 1986 Adam Fritz, 133 Main St., Afton, NY 13730 }
  2.  
  3. procedure vSGEFA ( var fID ;
  4.                     lda, n : integer ;
  5.                   var IPvt : IVECTOR ;
  6.               var InfoCode : integer ) ;
  7.  
  8.      { PROGRAM:                                            }
  9.      {                                                     }
  10.      {    SGEFA - Factors Real Single Precision Matrix     }
  11.      {            Using Gaussian Elimination.              }
  12.      {                                                     }
  13.      { VERSION:                 DATE:                      }
  14.      {                                                     }
  15.      {     1.0/TURBO Pascal 3.0    08/02/86                }
  16.      {     FORTRAN                 08/14/78                }
  17.      {                                                     }
  18.      { DESCRIPTION:                                        }
  19.      {                                                     }
  20.      {     SGEFA is usually called by SGECO, but can be    }
  21.      {     called directly with a saving in time if RCOND  }
  22.      {     is not needed.                                  }
  23.      {                                                     }
  24.      {     (time for SGECO) = (1 + 9/n)*(time for SGEFA)   }
  25.      {                                                     }
  26.      {     On entry;                                       }
  27.      {                                                     }
  28.      {        A     - Matrix to be factored.               }
  29.      {        lda   - Leading dimension of matrix A.       }
  30.      {        n     - Order of matrix A.                   }
  31.      {                                                     }
  32.      {     On exit;                                        }
  33.      {                                                     }
  34.      {        A     - Upper triangular matrix and multi-   }
  35.      {                pliers used to obtain it. Factoriza- }
  36.      {                tion can be written;                 }
  37.      {                   A = L * U                         }
  38.      {                where L is a product of permutation  }
  39.      {                and unit lower triangular matrices   }
  40.      {                and U is upper triangular.           }
  41.      {                                                     }
  42.      {        IPvt  - Pivot index vector.                  }
  43.      {                                                     }
  44.      {        InfoCode -                              }
  45.      {            = 0 Normal value.                        }
  46.      {            = k If  u[k,k] = 0.0. This is not an     }
  47.      {                error condition for this procedure,  }
  48.      {                but indicates that SGESL or SGEDI    }
  49.      {                will divide by zero if called.  Use  }
  50.      {                RCOND from SGECO for reliable indi-  }
  51.      {                cation of singularity.               }
  52.      {                                                     }
  53.      { SUBPROGRAMS:                                        }
  54.      {                                                     }
  55.      {    SAXPY       from BLAS                            }
  56.      {    SSCAL       from BLAS                            }
  57.      {    ISAMAX      from BLAS                            }
  58.      {                                                     }
  59.      { AUTHORS:                                            }
  60.      {                                                     }
  61.      {    Cleve Moler                                      }
  62.      {    University of New Mexico                         }
  63.      {    Argonne National Laboratories                    }
  64.      {                                                     }
  65.      {    PASCAL translation;                              }
  66.      {                                                     }
  67.      {    Adam Fritz                                       }
  68.      {    133 Main Street                                  }
  69.      {    Afton, New York 13730                            }
  70.  
  71. var
  72.      gID            : vARRAY Absolute fID ;
  73.      j, k, l        : integer ;
  74.      kp1, nm1       : integer ;
  75.      t              : real ;
  76.  
  77. begin
  78.                                 { Gaussian Elimination with Partial Pivoting }
  79.    InfoCode := 0 ;
  80.    if n > 0 then begin
  81.       nm1 := n - 1 ;
  82.       for k := 1 to nm1 do begin
  83.          kp1 := k + 1 ;
  84.                                 { Find l = Pivot Index }
  85.          iAk := VectorRead(gID,n,k,k,n-k+1,Ak) ;
  86.          l := ISAMAX(n-k+1, Ak[iAk], 1) ;
  87.          IPvt[k] := l + k - 1 ;
  88.                                 { Zero Pivot Implies Column Triangularized }
  89.          if Ak[iAk+l-1] <> 0.0 then begin
  90.                                 { Interchange if Necessary }
  91.             if l <> 1 then begin
  92.                t := Ak[iAk+l-1] ;
  93.                Ak[iAk+l-1] := Ak[iAk] ;
  94.                Ak[iAk] := t
  95.             end ;
  96.                                 { Compute Multipliers }
  97.             t := -1.0/Ak[iAk] ;
  98.             SSCAL(n-k, t, Ak[iAk+1], 1) ;
  99.             VectorWrite(gID,n,k,k,n-k+1,Ak) ;
  100.                                 { Row Elimination with Column Indexing }
  101.             for j := kp1 to n do begin
  102.                iAj := VectorRead(gID,n,k,j,n-k+1,Aj) ;
  103.                t := Aj[iAj+l-1] ;
  104.                if l <> 1 then begin
  105.                   Aj[iAj+l-1] := Aj[iAj] ;
  106.                   Aj[iAj] := t
  107.                end ;
  108.                SAXPY(n-k, t, Ak[iAk+1], 1, Aj[iAj+1], 1) ;
  109.                VectorWrite(gID,n,k,j,n-k+1,Aj)
  110.             end
  111.          end
  112.          else
  113.             InfoCode := k
  114.       end
  115.    end ;
  116.    IPvt[n] := n ;
  117.    if Aj[iAj+1] = 0.0 then
  118.       InfoCode := n
  119. end ;
  120.  
  121. { Copyright (C) 1986 Adam Fritz, 133 Main St., Afton, NY 13730 }
  122.